Universal properties of correlation transfer in integrate-and-fire neurons 
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One of the fundamental characteristics of a nonlinear system is how it transfers correlations in its 
inputs to correlations in its outputs. This is particularly important in the nervous system, where 
correlations between spiking neurons are prominent. Using linear response and asymptotic methods 
for pairs of unconnected integrate-and-fire (IF) neurons receiving white noise inputs, we show that 
this correlation transfer depends on the output spike firing rate in a strong, stereotyped manner, 
and is, surprisingly, almost independent of the interspike variance. For cells receiving heterogeneous 
inputs, we further show that correlation increases with the geometric mean spiking rate in the same 
stereotyped manner, greatly extending the generality of this relationship. We present an immediate 
consequence of this relationship for population coding via tuning curves. 
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Systems with spatially correlated, stochastic forcing 
provide a rich field of study, with applications ranging 
from polymer physics to population biology [lj. In neu- 
roscience, much theoretical work focusses on correlation 
(or synchrony) between coupled neurons receiving uncor- 
rected external fluctuations 2]. Nevertheless, in sensory 
systems, external fluctuating signals common to a group 
of unconnected cells 0, [l5| are a strong source of correla- 
tions. Moreover, signal-independent correlations can de- 
velop within strictly feedforward networks, due to conver- 
gent connections between layers [i, 14 1. Furthermore, the 
connectivity within biological neural networks is often 
sparse [H, suggesting that pairs of neurons rarely influ- 
ence one another directly [l2|. These observations have 
prompted experimental and theoretical efforts 

UMM to 

characterize the transfer of correlated input currents to 
correlated spikes for pairs of uncoupled neurons. Despite 
the straightforward setting of this problem, few common 
principles have emerged. 

Recently, we have observed that in response to fluctu- 
ating input currents with fixed levels of correlation, the 
spike (output) correlation between pairs of neurons in- 
creases with the spiking rate and was approximately in- 
dependent of the spiking variance (l3j . We provided em- 
pirical evidence for this correlation-rate relationship in in 
vitro cortical data and in IF models with 'hard' thresh- 
olds, for homogeneous cell pairs firing at low to medium 
rates. In this letter, we develop an analytical approach 
that both explains and greatly extends these observa- 
tions. Specifically, we obtain closed-form expressions re- 
lating output correlations to firing (spiking) rate, show 
that a refractory period implies the existence of optimal 
firing rates for correlation transfer, and show that corre- 
lation transfer depends approximately only on the geo- 
metric mean of the output firing rates for heterogeneous 
cell pairs. We demonstrate extensions to the quadratic 
integrate and fire model and to a broader analysis of 
the in vitro data. These results enable us to develop 
the first consequences of the correlation-rate relationship 



for cortical coding, where populations of heterogeneously 
tuned neurons relay information about the sensory envi- 
ronment. 

Model - The model consists of two IF neurons receiving 
common £(t) and independent £i(t) white noise inputs 
(Fig. UK): m 



rV' = f(Vi) + Mi + o-isTr [vT^&(t) + v^e(*)] , (1) 

i = 1,2. Here, Vi denotes the membrane voltage of cell 
i, t is the membrane time constant, and /(•) defines the 
subthreshold dynamics. We mostly take f(Vi) = —Vi 
(leaky IF), but in one section use f(Vi) = (quadratic 
IF). We adopt the standard threshold-spike-reset condi- 
tion: Vi{t + ) = Vr whenever Vi(t) = Vr; at such times, a 
spike is emitted, after which the voltage is held at Vr dur- 
ing the refractory period r r . Throughout, we set r = 1 
(we report rates in units of t _1 ). Of interest are the nor- 
malized output spike trains y%(t) — J^i^it ~ *f )> where 
t\ is the k th spike time of the i th neuron. The firing rate 
of the i th cell is V{ — (yi(t)}. 

The input current Ii(t) to cell i has a constant (DC) 
component p,i as well a stochastic component with ampli- 
tude o~i (Fig. UK) . The stochastic component to each cell 
is the sum of two independent gaussian processes, £,(t) 
andtft), with (&(t)&(i')) = 5 l3 8{t-t') and (&(t)£(f)> = 
= 1,2); note that £(t) is common to both cells. 
The weights \/l — c and yfc result in a correlation co- 
efficient c between Ii(t) and hit). The spike correla- 
tion coefficient, p, is computed @, H, EH EH by inte- 
grating the spike auto- and cross-correlation functions 
(yi(t)yj(t + T)) - ViVf. 
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Linear response calculations - Following methods 
of [IH, we derive an analytic expression for p valid 
to linear order in c; simulations of Eq. ([1]) verify that 



these methods produce good approximations for c up to 
~ 0.3 [13|, and were repeated to the same effect for the 
cases studied here (data not shown). 

The cross spectrum of the spike trains y\ (t) and y 2 (t) 
is P\ 2 (i-o) = {y*(v)y2 where the Fourier transform 
y i = -i= f Q T e lut {yi{t)~Vi)dt. Treating the common term 

y/c^(t)as a perturbation, we use the linear approxima- 
tion 



wM = y i>w , CTi VT^(w) + ^VcA Mj)(7 . VT ^(a;)|(a;) + 0(c). 

(3) 

Here „ f CT4a /t^(w) is the Fourier transform of a spike 
train obtained from a neuron driven only by the "private" 
input /ij+(Tj-\/i — c£i(t) (i.e., in the absence of a common 
stochastic input), and A^. 71-37(0;) is the susceptibility 



function of such a neuron [23 

Inserting Eq. ([3]) into the definition of the cross spec- 
trum we obtain: 

P 12 (lo) = C( r 1( 7 2 ^ i!(Ti (c)A AI2 ^ 2 (a;)P c (o;)+0( C 3 / 2 l4) 



where 



P ( (w) 



' i.f-ti ,(Ti — c ^j^j i&jV i.[±i .ai v^T- 

1,2. We also used the fact that the 0(c) terms in 
Eq. ([3]) arise from £, so are statistically independent of 
the 7^ CTiv 7 ! 3 ^' Moreover, we evaluate the susceptibility 
terms at c = 0, which is valid to the order of error shown. 
The Wiener-Khinchine Theorem yields Pi2(0) = 
Ci2(t)dt; also, the susceptibility to DC input is 
A WjO . 4 (0) = Thus, substituting Eq. g]) into Eq. © 
yields to lowest order in c: 



1. 



Here we used 

:f) = 0, i,j = 



ca±a 2 



dv\ du2 
dfj,! dfj.2 



CViCV 2 ^U^ 1 



(5) 



where we have additionally used the result for renewal 
point processes J Cu(t)<1t — CV 2 Vi] here, CVj is the 
coefficient of variation of the inter-spike interval distri- 
bution 0. 

Dependence of correlations on firing rate. - We first 
assume that the statistics of inputs to cells i and j are 
identical (/ij = /ij = /i and Cj = <7j = a), giving 



CV 2 l 



(6) 



Here S is the correlation susceptibility [l3j between 
the input current correlation c and the spike (output) 
correlation p. 

We explain three key features of Eq. ([6]) : i) S increases 
with firing rate v for low to moderate rates, ii) S is 
approximately a function of v alone, independent of the 
precise values of /1 and a or the CV. At low rates this 
relationship becomes exact. Hi) For large rates: S tends 
to a constant near 1 if r r = and tends to zero if r r > 0. 




FIG. 1: a) Schematic of two cells receiving heterogeneous, 
correlated white noise inputs, b) Rays with different /3 = fi/a; 
greyscale (color online) applies to the following panels, c) S 
vs. v increasing along the different rays for a limited range 
of rates for LIF cells, r r = 0. Here and for all figures with 
LIF model: Vt = 1, Vr = 0. Kramer's approximation shown 
as dotted line. Inset, on log-log axes, d) As (c), for wider 
range of rates. Dotted lines show Kramer's approximation at 
low rates and large fj,, a approximations at high rates. Sets 
of curves, from top to bottom: r r = 0, r r = 0.25, r r = 0.5. d) 
Similar, but for QIF model with r r = 0, Vt = 10, Vr = -10. 



To study the low rate case we treat the voltage of each 
cell as a particle in the potential 4>(V) — (V— p,) 2 /2. Low 
rates are obtained only when a < Vt and < a <C 1. 
In this limit, Kramer's theory can be applied to calculate 
the terms in S, by computing the mean (m) and the vari- 
ance (var) of hitting times of the threshold voltage Vr for 
particles starting at the bottom of the wells, and by not- 
ing that v = (m + r r ) _1 and CV 2 — z^ 2 var. Asymptotic 
results (Til give, to lowest order in a, 



and 



:exp(-a 2 ) = g(a) 



Z (2a -h) , where a = 



(7) 



Moreover, 



CV 2 = 1 to lowest order in a, as the threshold crossing 
times are exponentially distributed in this limit. There- 
fore, here we have S ~ v (2a — I /a) 2 . Eq. ([7]) may be 
uniquely inverted for a range of sufficiently small a (large 
a) to yield a — g (y). For small v we obtain: 



1 



(8) 
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so that asymptotically S is a function of v alone (Fig.QJ:). 

Asymptotics for large firing rate v are obtained from 
exact expressions for the firing rate and variance [20| . 
Using asymptotic approximations, we obtain for fixed a 
and large p 



S(fj,,a) 



(Vt - Vr) 

fJ,T r + (V T - V R ) 



Thus, the correlation susceptibility S approaches in the 
large p limit, as long as r r > 0. In the case r r = 0, we 
obtain S — > 1 as p — > oo (Fig. QJi). 

More generally, if p and a diverge jointly so that 
p/a = is constant (i.e., along a ray in Fig. \Tjp), we 
obtain S(p = (3a, a) ps Jg^Sv Here, Ki{f3) = 
(-2/V^+2/3e' 32 erfc( / 3)y _ „ 2 

[ Jf~ Sertfvl and K *® = ^ ^ US " 
ing asymptotic expressions for v at large p and cr, we 
obtain 



S(fl, cr) PS = ^T^j (1 - 



(9) 



The constant Ki(/3) / K2(f3) increases from 0.918 to 1, as 
(3 increases from to oo. Therefore, if t t = 0, S(y) is 
asymptotic to values in the interval [0.918, 1] for large v, 
as shown in Fig. [TJi. Furthermore, if t t ^ then S(y) 
tends to zero at v — r,? 1 , with slope in [—1, —0.918]. 

In summary, we have shown that, at low rates, S is 
asymptotically a function of v only, regardless of specific 
values of a and p and hence of the CV. At high rates, 
if r r = then S is bounded in a neighborhood of 1; if 
r r > then S tends to zero within a narrow "cone" at 
rates near the limit v — t^ 1 . It follows that a refractory 
period implies that S obtains a maximum at a partic- 
ular firing rate, at which the neuron is "best-tuned" to 
transfer input correlations to output correlations. For in- 
termediate rates, numerical evaluations of Eq. © show 
that S remains an approximate function of v alone, a fact 
that greatly simplifies the statistical description of the 
LIF model. We note that these findings hold more gen- 
erally: FigQJ shows analogous results for the quadratic 
integrate and fire model, demonstrating that our results 
depend on neither the linear subthreshold dynamics nor 
the "hard" threshold condition of the LIF model. 

Effects of heterogeneity. - We return to the general 
Eq. (0. In neuroscience, several studies plot the out- 
put correlation of two cells versus the geometric mean of 
their firing rates ^/uiU2 [1,0, E3]- Fig. [2^ shows the corre- 
sponding predictions of Eq. ([5]), as the input parameters 
{pi, o\, p2, 02} are sampled over a wide range, producing 
very different combinations of v%, CVi, 2 = 1,2 (here and 
below we take r r = 0). Remarkably, when p is plotted 
vs. y/v\vi, we see approximately the same correlation- 
rate relationship as for the homogeneous case. At low 
rates, this relationship becomes almost independent of 



the ratio h — vxjv^ (y\ > V2) between the cells' rates; at 
higher rates, it depends only weakly on h. Thus, in the 
limit of small c we have: 



p = cy S(/ii,o-i)S(fJ,j,o-j) « c\J S{vi)S{vj) PS cS{y/ViVj). 

(10) 

Here, the first approximation follows from our results for 
the homogenous case, and the second follows empirically 
from Fig. [2^l. 

To see how surprising this is, consider the four volt- 
age traces in Fig. [5^, each of which displays either a 
high or low CV, and a high or low rate (y =0.47 or 
0.047, resp.). Nevertheless, any combination of these 
cases that gives the same geometric mean (AB, AC, BD, 
CD, y/v\V2 = 0.15) gives the same value of S = 0.55, 
within less than 1% (square on graph). This is despite 
the strongly distinct structure of the correlation func- 
tions for these different pairings (AB, CD shown). Other 
striking predictions similarly follow from Eq. (fTt)|) : for 
example, simply adding a DC current to one of a pair of 
cells (thus increasing, e.g., vC) should increase their spike 
correlation. 

Why does Eq. ^ hold? Fig.[T]shows that S{n u a t ) ps 
S^), with exact agreement at low rates. The approxi- 
mation J S{vi)S{vj) ps S(«yUiUj) would be exact if S(ui) 

was a power law, S{vi) cx vj , and Fig. [Tt (inset) shows 
that this is an excellent approximation at low rates. 
Hence Eq. (fT0|) holds precisely at low rates, so that the 
different areas in Fig. [5^ converge. To explore when the 

approximation J S{vi)S{vj) ps S{^JT^vj) breaks down, 
we study the special case where the ratio /3j = Cj j pi — 1 
is identical for both cells (Fig. [2Jd). Fig. [2J: shows level 
sets of this surface on log-log axes (thick lines), which al- 
most superpose with level sets of yJvpTj (thin black lines) 
for low rates, or for values of h reasonably close to 1 
(shaded areas). For larger rates, or greater heterogene- 
ity, these level sets start to diverge. However, even then 
the disagreement is quite mild: for h up to 10, values 



of S along level sets of 



do not vary by more than 



ps 10%, because S varies only weakly with large rates. 
This agreement stands in contrast to that between level 
sets of S and of the arithmetic mean (i>\ + 1/2) /2, shown 
as grey lines (red online), indicating that the choice of 
the geometric mean is somewhat special in capturing the 
rate-dependence of S. Fig. [JJi shows that this relation- 
ship holds not only for the LIF model, but also for our in 
vitro experiments. Here, pairs of cortical cells in mouse 
brain slices were injected with correlated, noisy, low-pass 
filtered currents, replicating the configuration of Fig. la 
(see [3| for experimental methods). Output spikes were 
recorded via whole-cell patch clamp, and p was computed 
from paired spike trains as described above. 

For an illustration of the functional consequences of the 
rate-correlation relationship for heterogeneous cells, con- 
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FIG. 2: (a) S vs. ^Jv\ V2, from Eq. Input parameters 

{/ii , tTi, /i2, o"2} were sampled over a 4-D cube (fiij £ [.2,8.2], 
Ci,j G [.2,8.2], grid spacing 0.01-.05, with denser grid for 
/i, a < 1.6,2.2). Overlapping regions shaded from dark to 
light contain all points with levels of h = u\/v2 varying over 
increasingly broad ranges. Square marks value for the pairings 
AB,AC,BD,CD shown in the inset (voltage traces and auto- 
and cross-correlation functions plotted); see text. Horizontal 
scale bar: 2r for voltage traces, 4t for correlation functions. 
Vertical scale bar: .25 for correlation functions, .5 for traces 
B,D, 1 for traces A,C. (b) Surface/contour plot of S, where 
a — n for both cells, (c) Contour plot on log-log axes; see 
text. (d,e) S vs. ^fv\V2 for in vitro cortical experiments; (d), 
raw data 1269 pairwise comparisons obtained from 20 cells 
using c = 0.5, (e), population average and standard deviation 
to show trends. 
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FIG. 3: (a) Firing rate in response to a given stimulus for a 
population of neurons with tuning curves evenly distributed 
across the stimulus (8) space, (b) The resulting crosscorrela- 
tions pij for this stimulus, from Eq. ((5]). 



cal [13] studies, follows directly from our results, without 
requiring a separate assumption that the input correla- 
tion is itself tuned. We also predict stronger correla- 
tions for neurons with "preferred" tuning to a stimulus, 
as in llfJl . Such an extension of stimulus selectivity from 
rate [2l| to correlation tuning has broad and novel im- 
plications for information encoding in spiking neurons. 
Further, correlations in output spike trains are effective 
in propagating activity downstream [4| , thereby facilitat- 
ing the decoding of correlation-based information. 

Our methods can be extended: Eq. |(SJ| holds for any 
cell or cell model producing spikes as a renewal process, 
and the preceding forms of this expression are valid even 



for nonrenewal spiking (enabling, e.g., analysis of tempo- 
rally correlated inputs pj). Moreover, our formulas and 
their reduction to a compact correlation-rate relation- 
ship may be a critical component for mean field models 
of sparsely coupled neural populations, where treatment 
of second-order statistics has previously been restricted 
to spike train variance (e.g., [221]). A mean field frame- 
work that rigorously includes correlations would connect 
neural population dynamics @, 0, [22| with the emerging 
role for correlations in sensory coding. 
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sider a population of neurons whose rates Vi are tuned to 
a sensory variable 6 as in Fig. 17|, l2l| . The correlation 
Pij between pairs of these cells approximately depends on 
their geometric mean rate yjv\vi according to Eq. (|10p . 
producing the correlation tuning shown in Fig.[3jD (input 
correlations are taken to be constant: cy = c). Thus, the 
fact that only neurons with similar tuning display sub- 
stantial correlation, as in empirical 1^, II] and theoreti- 
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